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TWO-DIMENSIONAL ANALYSIS OF HEAT AND 

MASS TRANSFER IN POROUS MEDIA USING THE 

STRONGLY IMPLICIT PROCEDURE 

By Donald M. Curry 
Lyndon B. Johnson Space Center 

SUMMARY 


The conservation relations describing the transient heat and mass transfer in a 
porous matrix are coupled, nonlinear partial differential equations. The equations are 
solved in finite difference form by using a recently developed numerical technique (the 
strongly implicit procedure). The application of the procedure is discussed. Sample 
numerical results demonstrate the two-dimensional characteristics of heat and mass 
transfer in a porous matrix. By imposing external pressure distributions, external 
heat flux distributions, and backface mass flux injection distributions, the effects on 
the internal pressure distributions, the internal mass flow rates, and the internal tem- 
perature distributions are obtained. 


INTRODUCTION 


An understanding of heat transfer and fluid flow in porous media is important in 
many engineering fields, such as soil mechanics, powder metallurgy, chemical proc- 
essing, and petroleum reservoir recovery. The successful application of ablation ma- 
terials or transpiration cooling systems to protect ballistic missiles and manned 
spacecraft has motivated investigations of the heat- and mass-transfer processes oc- 
curring in a porous matrix at high temperatures. The basic thermal features of an 
ablative system are a high surface temperature, a steep temperature gradient within 
the material, and the transpiration of gases (generated by thermal decomposition) 
through a high-temperature porous matrix. The high surface temperature permits ap- 
proximately 80 percent of the external heat flux to be reradiated by the surface. The 
steep temperature gradient permits a cool backface with a hot external surface. The 
transpiring gases remove some heat being conducted into the material and also thicken 
the boundary layer, thereby decreasing the convective heating to the wall. 

Both ground-based and flight test data (obtained from a charring ablator used on 
the Apollo spacecraft) have indicated a need to extend the currently available analytical 
techniques to describe more completely and to predict the expected performance of the 
thermal protection system. An analysis of the thermal response of the Apollo ablative 
heat shield (ref. 1) has shown a discrepancy in current design methods for calculating 



surface recession. The prediction of the one-dimensional analytical model (ref. 1) 
compares favorably with flight data but underestimates the surface recession results 
in ground-based tests. Improvement of the prediction technique requires further study 
of phenomena associated with aerothermodynamic environment, material chemistry, 
and multidimensional heat and mass transfer in a porous char matrix. This report con- 
siders only two-dimensional heat and mass transfer in a porous material. An analysis 
by Bush and Dow (ref. 2) of the pressure fields within cylindrical and hemispherical 
ablative chars under steady-state conditions has indicated that the gas flow through the 
char can be significantly different from that predicted on the basis of a one-dimensional 
model. The theoretical results further indicated the possibility of an influx of boundary- 
layer gas into the char. However, Bush and Dow considered no coupling of the energy 
and momentum equations in the analysis. 

Heretofore, analytical and experimental studies in the literature are generally 
not concerned with the coupled heat-transfer/fluid-flow problems that occur in the chem- 
ical, petroleum, and aerospace industries . However, Bland (ref. 3) showed the exist- 
ence of coupling between the momentum and energy relationships when the steady-state, 
Darcy continuity equation was used to describe the flow of a gas through a porous 
medium. Bland considered a temperature difference between the gas and solid and heat 
conduction in the solid but neglected the gas heat conduction and storage in the fluid 
energy equation. The coupled equations are developed for the case of three independent 
space variables and one time variable. The theoretical treatment of Bland shows clearly 
the coupling of the momentum and energy equations used to describe simultaneous heat 
transfer and fluid flow in a porous matrix. 

More recently, Schneider et al. (ref. 4) performed a two-dimensional, steady- 
state analysis of an active transpiration cooling system. For an incompressible coolant 
flow through a porous nose tip, the two-dimensional solutipn was found to be important 
in assessing the internal coolant distribution and in determining the driving pressure 
required to prevent coolant starvation at the nose tip. For the purpose of comparison, 
Schneider et al. performed a one-dimensional analysis that resulted in a nonconservative 
design with respect to expulsion system pressure requirements. 

Del Casal (ref. 5) investigated the effects of multidimensional flow on heat trans- 
fer for the stagnation region of blunt axisymmetric surfaces. The internal-flow solu- 
tions were coupled with those of the external boundary layer by matching solutions at 
the outer surface of the porous matrix. Although only steady- state solutions were per- 
formed, Del Casal included the effects of inertia in the momentum equation; these had 
been neglected by Schneider et al. 

Pittman and Howser (ref. 6) analyzed the transient response of an axisymmetric 
charring ablator, considering internal flow. They assumed thermal equilibrium between 
the gas and solid phases. Surface recession and pyrolysis interface rates were shown 
to be affected by multidimensional flow of the pyrolysis gases through the char layer. 

Numerous analytical models for predicting pressure distribution, mass flow rate, 
and thermal response of porous materials are published for a wide variety of boundary 
conditions. One predominant assumption in a majority of these analyses is the inde- 
pendence of the energy and momentum conservation relations used for performance 
predictions. Voluminous material is available on both heat transfer and pressure-drop 
flow studies for a porous matrix. However, the use of porous materials in spacecraft 
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thermal protection systems has emphasized the importance of understanding the more 
complex, coupled relationship. The general lack of available solutions that consider 
the coupled energy and momentum relations used in describing the thermal response 
of a transpiration cooling system motivated this study. 

SYMBOLS 


A,B,C,D,E,G,H,Q 

AX, AY 

b,c,d,e,f 

C 


g 

h 

v 

K 

k 


k ,k 
x’ y 

L , L 
x’ y 


parameters known from previous time level and previous 
iteration 

functions defined in equation (B3) 
coefficients defined in equation (21) 
specific heat 

gravitational constant 
volumetric heat-transfer coefficient 

proportionality factor 
thermal conductivity 

thermal conductivity in x and y directions, respectively 
length in x and y directions, respectively 


rh 

m m 

npx 

npy 

P 


mass flow rate, pV 

volume source term 

total number of nodes in x direction 

total number of nodes in y direction 

pressure 


q convective heating rate 

q’" volumetric heat source (or sink) 

R universal gas constant 

T temperature 
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radiation equilibrium temperature 


eq 


1 

T » 

t 


initial temperature 

u nkn own temperature in difference equations 
time 


U total internal energy 

u velocity in x direction 

V superficial velocity 

Vj j intermediate coefficient used in the solution algorithm 

v velocity in y direction 


v 


x,y 

a 

0 

r 

5 


6 

M 

P 

a 

T 

$ 

* 


velocity in a pore 
space coordinates 
viscous coefficient 
inertial coefficient 
iteration parameter 

dependent-variable difference between the new time step or 
iteration and the previous time step or iteration 

emissivity 

viscosity 

density 

Stefan- Boltzmann constant 
viscous stress tensor 
viscous dissipation 

porosity, ratio of void area to total area 
convective heat-blocking ratio due to mass injection 
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Subscripts : 


e 

f 

i 

i 

o 

Ref 

Stag 

s 

w 


external 

fluid 

x-direction node location 

y-direction node location 

without blockage 

reference 

stagnation 

solid 

matrix surface 


Superscript: 

n iteration or time step 


Operators : 

(') indicates parameter at next time step 

v differential operator 


THEORETICAL RELATIONS 


The unsteady heat transfer of a compressible gas flowing within a porous matrix 
can be described by the solution of the conservation equations for mass, momentum, 
and energy. The fundamental assumption is the use of a continuum model for the porous 
matrix. Briefly, the volume element used in deriving the conservation equations is 
assumed large with respect to the individual pores so that the flow field may be con- 
sidered a continuum. The application of the conservation equations for mass, momen- 
tum, and energy to the unsteady flow of a gas, with heat transfer, through a porous 
matrix is presented in appendix A. 
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For the case of a two-dimensional heat transfer and fluid flow in a porous medium, 
the governing differential relations can be written as follows. The Darcy continuity 
equation is 



where $ = porosity 
P = pressure 
t = time 

T = temperature 

x,y = space coordinates 
a = viscous coefficient 
M = viscosity 
ft = inertial coefficient 
p = density 

u = velocity in x direction 
v = velocity in y direction 
R = universal gas constant 
m ,M = volume source term 

The energy considerations when Tq u ^ is not equal to T ^ are as follows. The 
fluid energy equation is 



(2) 
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where = specific heat 

k = thermal conductivity 

h = volumetric heat-transfer coefficient 
v 

q m = volumetric heat source (or sink) 

The solid energy equation is 


ST a f 

pc -^ = ~k 

s p 3t Sx ' 
*s 



(1 


^( T s- 


T f ) + q’ 


(3) 


These three partial differential equations plus two algebraic relations, the Darcy 
pressure drop (eq. (All)) and the perfect gas equation of state (eq. (A12)), form the 
complete set of equations. 


BOUNDARY AND INITIAL CONDITIONS 


The following specification of boundary and initial conditions is an important phase 
in describing the physical system being modeled. In the present analyses, the physical 
domain consists of a rectangular matrix with dimensions L x and as shown in 

figure 1. When these analyses are applied to the case of a charring ablator, the ablator 
exterior surface is represented by the plane x = 0 and the char thickness is L x - The 

fluid enters the matrix at the rear face x = L . The mass flow rate at the center rear 

X 

face m(L x ,0) is used as the base value to nondimensionalize local mass flow rates. 

All injected gases flow toward the heated surface x = 0; no flow reversal is permitted. 
Mathematically, the boundary conditions for the matrix and fluid energy equation at the 
exterior surface are 



Figure 1.- Two-dimensional porous 
matrix. 


x = 0; 


ST. 


-k 


s 


s Sx 


^net 



S 2 T, 


Sx 


= 0 


(4) 

(5) 
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where q = convective heating rate 

a = Stefan- Boltzmann constant 
e = emissivity 

The net heat at the exterior surface is specified by equation (4). Here, p is defined 
as the ratio of the convective heat flux with blowing to the convective heat flux that would 
exist without mass injection. For the boundary condition as shown by equation (5), the 
matrix surface does not require that fluid and matrix temperatures be equal; it does 
not preclude T = T f . This boundary condition permits conduction in the fluid at the 
fluid exit. s 

The temperatures of the gas and solid are assumed to be equal along the injection 
plane. This is consistent with the gas being formed at a local solid decomposition 
temperature. 


X = V T f = T s 


( 6 ) 


A thermally adiabatic back surface exists along the plane of mass injection. 


9T 

x = L x ; 


= 0 


(7) 


Symmetric boundaries are assumed in the y direction. 


y = 0 ’V 


ST f 

~w 


= 0 , 


3T c 

~9y 


= 0 


( 8 ) 


The unsteady boundary conditions on the momentum equation permit a variable pressure 
in the y direction. 


x = 0; P(0,y,t) = P e (y,t) 


(9) 


x = L ; = 0 

x 3x 


(10) 
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where P g is external pressure. The boundary condition illustrated by equation (10) 

does not permit flow reversal at the injection plane. This is a reasonable assumption 
since an impermeable seal exists a small distance from the injection zone. Along the 
y boundary 


y = 0; 


9P 

9y 


= o 


y = L ; 

* y 


§ = 0 

9y 


(11) 

( 12 ) 


The initial pressure distribution is assumed constant and equal to the local external 
pressure at the start of the transient solution. No flow can occur in the y direction 
along the x = L boundary, with symmetrical flow assumed at the boundaries for 

X 

y = o, y = L y . 

Equations (1), (2), and (3) describe the flow and heat transfer in a porous matrix 
and are coupled, nonlinear partial differential equations. The solution is achieved by 
approximating the partial derivatives by use of suitable finite difference expressions 
involving the independent and dependent variables. This procedure leads to a set of 
algebraic equations in the dependent variable from which the final value can be deter- 
mined. Three basic finite difference formulas (forward, backward, and central differ- 
ence) can be used to obtain the necessary relations. Once the finite differencing method 
has been chosen, a variety of solution techniques is available, each of which has its own 
peculiar advantages and limitations. For example, an explicit technique may require 
small integration time steps to ensure stability; the implicit technique seemingly does 
not have a time-step limitation but may require the solution of a large algebraic matrix. 

In this investigation, an implicit central- difference numerical technique is used 
to obtain solutions for the momentum, energy, and continuity relations; it requires 
evaluation of the coefficients of the dependent variable at the previous time level and/or 
iteration as well as the solution of a set of linear difference equations for the dependent 
variables. 


IMPLICIT DIFFERENCE 
EQUATION FORMULATION 


A general two-dimensional region 
with an x-y grid system imposed is shown 
in figure 2. The implicit difference for- 
mulation used to solve equations (1), (2), 
and (3) can be illustrated conceptually 



Figure 2.- Two-dimensional region 
with rectangular grid. 
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by consideration of the simpler problem of pure conduction. The transient heat con- 
duction equation for this two-dimensional region is 



(13) 


An implicit finite difference equation for each grid point (i, j) within the 
specified region can be written as 



where, for the purposes of this illustration, k x and have been considered constant. 

The prime on each temperature represents an unknown temperature. Equation (14) can 
be written as 


i,3 


(Ay Y 




k 

X. 


1,3 

T’ - 

.(ai?. 

1-1,3 

"k 
x. . 


x ?3 

T T 4- 

.(Ax) 2 . 

1+1,3 


, 2k 2k 

x. . y. . 
-P + — bj. + -114 

At , . \2 , A sZ 

(Ax) (Ay) 




PC 


T.' . 1 = - 
1,3+1 


T! . 


T. . 


At i,j 


(15) 


Rewriting equation (15) yields 


AT.’ . . + BT.' . + CT.' . + DT.’ . . + ET.’ . , = Q. . 
1,3-1 1-1,3 1,3 1+1,3 1,3+1 +,3 


(16) 


Equation (16) has five unknown temperatures per grid point (i, j). The values of 
A, B, C, D, E, and Q are known based on the previous time level and previous 
iteration. A set of equations similar to equation (16) can be written for all (i, j) grid 
points within the region and on the boundaries. This matrix of equations then can be 

inverted to yield the unknowns T.’ .. For large systems of equations, this matrix 

1,3 

solution can become very time consuming. 

The development of the alternating direct implicit (ADI) method for the two- 
dimensional system proposed by Peaceman and Rachford (ref. 7) reduced the number 
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of unknowns to three, as obtained for simple one-dimensional problems. Basically, the 
ADI solves the equations in one direction while assuming the dependent variable in the 
second dimension to be constant over the time interval. As an example, consider equa- 
tion (16) for the first time step, in the x direction. 


BT.' 

l- 


, . + CT.’ . + DT.’ . = Q, . - AT. . , - ET. . 1 
1,3 i,3 1+1,3 i,3 1,3-1 1,3+1 


(17) 


The temperatures T^ ^ and T^ ^ are known from the previous time step. Equa- 
tion (17) yields a tridiagonal matrix. 



+ C 1 T 2,3 


= q i,j 

Vi,) 

+ B 2 T 2,j 

+ C 2 T 3,j 

= Q 2,i 


A 3 T 2, j 

+ B 3 T 3,j +C 3 T 4,j 




A N-l T N-2, j +B N-l T N-l,j +C N-l T N,j 

“ Q N-l,j 



A N T N-1, j +B N. T N,3 

= Q N,j 


Equation (18) is solved for the unknown temperatures by Gaussian reduction. Then, for 
the second time step, writing equation (16) in the y direction yields 


AT.’ 


i,3 


, + CT.’ . + ET.’ . n = Q. . - BT. . . - DT. , . 
■1 l,] i,]+l ",3 1-1,3 1+1,3 


(19) 


Once again, a tridiagonal matrix as shown by equation (18) is obtained. The advantages 
of solving a tridiagonal matrix over the matrix generated by equation (16) are apparent. 
The ADI method is not used in this study; rather, the strongly implicit procedure (SIP) 
is used in the next section. 


STRONGLY IMPLICIT PROCEDURE 


Recently, Stone (ref. 8) developed a new iterative method for solving sets of alge- 
braic equations that occur for multidimensional systems. This method is known as the 
strongly implicit procedure and has been used successfully in solving the two-dimensional, 
steady -state heat conduction problem as well as multidimensional flow in a petroleum 
reservoir. The foundation of the SIP calculation method lies in the approximate factor- 
ing of the five-diagonal matrix (five nonzero elements in each row of matrix) generated 
by equation (16) into three-diagonal upper and lower triangular matrices. The detailed 
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mathematical matrix reduction process, required to derive the upper and lower trian- 
gular matrices, is presented by Stone (ref. 8) and more recently by Bozeman (ref. 9). 

The process of defining the upper and lower matrices involves the temperatures 
of the grid point (i, j). As shown in figure 3, the modified matrix equation involves not 
only the grid points indicated by the solid circles but also the temperatures at grid 
points (i- l,j+l) and (i+l,j-l) shown by the open circles. 

The influence of these new temperatures in the modified equation is balanced par- 
tially by subtracting approximately equal terms. The modified equation for the unknown 
side of equation (16) becomes 



H 


Figure 3. - Grid model for definition of 
upper and lower matrices. 


+ Cjt. 


a. .t: . , + b. .t* . + c. .t: . 

x,3 1,3-1 1,3 1-1,3 1,3 1,3 

+ d. .t.’ . . + e. .t: . . +h. .[t: , . . 

1,3 1+1,3 i,J 1,1+1 1,3 L 1+1, 1-1 

- rf-T! . + T! . . + T.' . . VI 

\ 1,1 1+1,1 1,1-1/J 

= Q l,i 

( 20 ) 


where y is the iteration parameter. The first five terms are simply the original 

left side of equation (16), and the H. G. . are used to cancel partially the temp- 

1, 1 i, 3 

eratures at (i-l,j+l) and (i+l,j-l). 

The equations used in this algorithm for SIP are 


A. . 

b. . 

1,1 l+re M _j 


c. . = 


B. . 

1,1 


%j 1 + 


H. . = b. .e. . . 
i,3 i,l i,3-l 


G. . = c. .1 . . 
i,l i,3 i-l,3 


(21a) 

(21b) 
(21c) 
(2 Id) 
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(21e) 


d i J = C i,i + rH i,i + y °i,r \ j 


D. . - yH. . 
P - hi i,3 
i.j d. . 

f E i,j ~ yG j,j 

^ “ d. }j 


(21f ) 


(21g) 


and 


R ” i ■ S,J - M.i-1 + + c i^,J + D M T i + i,J + E i,J T i l ,i+i) 


( 22 ) 


where the exponent n is the iteration or time step. After solving for the coefficients 
b, c, d, e, and f, then all the intermediate coefficient values of . are solved by 


b. .V. . , + c. .V. . . + d. .V. . = r: . 
1,3 1,3-1 1,3 1-1,3 1,3 i,3 i,3 


(23) 


Equations (21), (22), and (23) consti- 
tute the algorithm to be solved for the un- 
known variable t!* + . . These equations are 

valid for all points within the grid system. 
The following conditions, which must be 
satisfied for these equations on the bound- 
aries, are illustrated in figure 4. 


E = 0 


B = 0 


0 = 0 


A = 0 


The coefficients generated by equa- 
tions (21), (22), and (23) are intermediate 
coefficients used during the elimination 

process. The V. . quantities are com- 
i,3 

puted from equations (21), (22), and (23) 
by starting at the line j = 1 and letting i 
take on values from 1, 2, 3 ... npx nodes. 
Then j is incremented by 1 and the proc- 
ess is repeated. 


x,i 


A. n ; E. =0 for y =0,y = L 
i,0 i.npy ’ ’ y 


B„ 0 . = 0 for x = 0,x * L v 

0,J npx,j x 


Figure 4.- Illustration of boundary 
condition restrictions . 
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After all 

solved for 6? + . 
j = npy. 1,3 


V. . coefficients have been generated and stored, equation (24) 
1,3 

by simply reversing the order; that is, starting with i = npx 


is 

and 


R n+1 TJ . „ R n+1 . R n+1 
6. . = V. . - e. .6. , . - f. .6. . 1 
i,3 1,3 i,3 1+1,3 i,3 1,3+1 


(24) 


, R n+1 _n+l m n 

where 5 = T - T 

n = iteration or time step 


Therefore, the values of temperature at the new time level (n+1) can be calculated 
from a knowledge of T at the current time level (n). Stone (ref. 8) recommends a 
variation of the procedure by carrying out the calculation just described in reverse 
order on each alternate time step. This has the effect of using the temperatures 
T. . . . and T. . . , on the odd steps. 

i-l, 3+1 i+l, 3-1 

The value of the iteration parameter y lies between 0 and 1 . The minimum value 
is not critical and could be zero; however, Stone states that the maximum value can be 
critical for obtaining convergence. For a heat conduction problem with constant prop- 
erties, Stone recommends 


y = 1 - 
max 


min 


2(Ax)' ! 


k (Ax)' 1 

1 +-* — . 


k x (Ay)' 


2(Ay) 2 

k x (4y) 2 

1 + o 

k y (Ax) 


(25) 


It should be pointed out that equation (25) is for a particular problem, and it may be 
necessary to vary y for other types of problems to determine its effect. In this work, 
a constant y value of 0. 95 was used successfully. The SIP formulation as applied to 
a simplified form of the fluid energy equation (2) is presented in appendix B. 


SOLUTION TECHNIQUE 


The two-dimensional nonlinear conservation equations (1), (2), and (3) were ap- 
proximated by a set of finite difference equations compatible with Stone’s algorithm. 
These equations are linearized by holding the nonlinear coefficients constant for each 
iterative calculation. The coefficients are then updated after each iteration. This 
set of linear equations is solved by application of equations (21) to (24). Because the 
equations are coupled through the nonlinear coefficients, the following iterative scheme 
was used as shown in figure 5. 
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Figure 5.- Numerical iteration 
procedure. 


The procedure is (1) to solve the solid 
energy equation by using the previous level 
values for flow rates, temperature, and 
physical properties; (2) to solve the Darcy 
continuity equation to obtain new pressures 
and velocities; (3) to obtain the fluid energy 
equation solution by using these latest values 
for flow rate, pressure, and temperature to 
obtain new fluid temperatures; and (4) to 
perform an independent calculation of the 
continuity equation (eq. (Al)) to ensure mass 
conservation. Three iterations per time 
step normally were required to obtain a 
satisfactory mass balance. The largest 


error occurs at the surface node where the pressure is being specified. Gas density 


for this node is evaluated from the continuity equation instead of the perfect gas law. 


NUMERICAL RESULTS 


The solution is capable of describing 
a wide variety of physical situations; how- 
ever, only a limited number of representa- 
tive cases are presented. More complete 
parametric comparisons for the porous ma- 
trix are provided in reference 10. Numeri- 
cal results are provided for cases similar 
to those experienced by the Apollo configu- 
ration. The external pressure distribution, 
the external heat flux distribution, and the 
backface mass flux injection distribution 
are shown in dimensionless form in fig- 
ure 6. The external distributions (mass 
flux injection at the rear boundary x = L , 

pressure and heat flux at the external bound- 
ary x = 0) were not time dependent in the 
cases under discussion. To identify the 
flow characteristics in the porous media, 
the decoupled solution (isothermal solution) 
is considered first. The effects of heat 
transfer on the internal flow characteristics 
are then considered by solving the coupled, 
two-dimensional equations. 



Figure 6. - External pressure gradient, 
heat flux, and mass injection flux 
ratios used to simulate typical Apollo 
conditions; where Pgt ag = p (0> 0), 

%ag = i!<0 ’ 0) ’ a ” d *Ref = lh(I V 0) - 
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Isothermal Solutions 


A metallic matrix material is considered for examples of isothermal solutions. 

For the metallic matrix, the porosity <p is 0.42, the viscous coefficient a is 

10 -9 4 -1 

1.619 x 10 m , and the inertial coefficient (3 is 1.634 x 10 m . 


If the Darcy continuity equation (eq. (1)) is solved for an isothermal case, the 
effect of the external pressure field on the internal flow in the matrix is demonstrated. 
An external pressure distribution (fig. 6) and a temperature of 556 K were used in the 
isothermal solution. A constant, uniform mass flux input at the rear surface was used. 
Air was used as the fluid flowing through the metallic matrix material and ablative char 
matrices; the corresponding properties of air were obtained from Kreith (ref. 11). 

3 2 

Based on an external stagnation pressure of 2. 03 x 10 N/m ', the local mass 
flow-rate distribution in the x direction is shown in figure 7. The mass flow rate has 


been normalized by dividing by the mass injection at the backface 


m = m(L ,0). 
Ref x’ 


Because the dimensionless mass flow rate does not vary substantially from unity, the 
internal flow is essentially one dimensional, indicating little effect of the external pres- 
sure distribution. The internal pressure distribution in the matrix is shown in figure 8. 



Figure 7.- Mass flux ratio normal to 
the surface for a two-dimensional, 
steady isothermal flow with the ex- 
ternal pressure gradient (metallic 
matrix) ; where Pg ta g = 2 . 03 x 

10^ N/m^ and ihj^ = HO. 3 kg/m^-hr. 


Figure 8.- Internal pressure distribu- 
tion for a two-dimensional, steady 
isothermal flow with the external 
pressure gradient (metallic matrix); 


where P c , = 2. 03 x 10 3 N/m 2 
Stag 2 

m Re f = 110.3 kg/m -hr. 


and 
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However, additional results with higher stagnation pressures show significant departure 
from one -dimensional flow within the porous matrix. For the case of an external stag- 

3 2 

nation pressure of 6. 08 x 10 N/m , substantial variation in the mass flow-rate distri- 
bution (fig. 9) is in the stagnation region (small y/L ). The corresponding internal 
pressure distribution is shown in figure 10. y 



Figure 9. - Mass flux ratio normal to 
the surface for a two-dimensional, 


steady isothermal flow with the ex- 
ternal pressure gradient (metallic 


matrix); where 
10 3 N/m 2 and 


P a , = 6.08 x 
Stag 2 

liiRef = 109. 6 kg/m -hr. 


Figure 10. - Internal pressure distribu- 
tion for a two-dimensional, steady 
isothermal flow with the external 
pressure gradient (metallic matrix); 


where Pg ta g = 6. 08 x 10 3 N/m 2 
^Ref = ® kg/m 2 -hr. 


and 


Nonisothermal Solutions 

For the nonisothermal case, the combined effects of a variable pressure and heat 
flux acting on the external surface (0,y) and of a variable mass injection rate at the 
backface (L ,y) are obtained by using the dimensionless distributions of figure 6. The 

coupled results shown in the following figures were obtained for the conditions of both 
thermal equilibrium and thermal nonequilibrium between the solid and fluid phases. In 
these nonisothermal considerations, a porous char matrix has been used. The thermo- 
physical properties of the char material are shown in table I. All heat-transfer calcu- 

2 

lations are for a stagnation point heating rate of 567. 5 kW/m and an initial temperature 
of 556 K. 
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The temperature distribution has been nondimensionalized by using the equilibrium 
surface temperature. The effect of the variable external heat flux on the temperature 
distribution demonstrates the two-dimensional effects as indicated in figures 11 to 14. 
The temperature distributions shown in figures 11 and 12 were obtained by using a vari- 
able surface heat flux distribution with a constant mass injection rate at the rear bound- 
ary (L ,y). The effect of the variable heat flux on the in-depth mass flux distribution 

(fig. 13) emphasizes the two-dimensional nature of the internal flow caused by in-depth 
pressure and density distributions. The internal pressure distribution is shown in 
figure 14. 


TABLE I.- THERMOPHYSIC AL PROPERTIES 
OF POROUS CHAR MATRIX 
(a) Constant property values 


Property 

Value 

3 

Density, p, kg/m 

Surface emissivity, £ 

Porosity, 0 

-2 

Viscous coefficient, a, m ... 
Inertial coefficient, /3, m' 1 ... 

320 
0.75 
0.66 
6.22 x 10 10 
6.31 x 10 5 


(b) Thermal conductivity k g as a function of temperature 


Temperature, K 

k s , W/m-K 

256 

0.00055 

422 

.00084 

811 

.00084 

1033 

.00175 

1144 

.00337 

1478 

.00597 

1589 

.00631 

2172 

.00914 

2200 

.00956 

2367 

.00956 

2700 

.00751 

2783 

.00685 


(c) Specific heat C as a function of temperature 
p s 


Temperature, K 

C , J/kg-K 
p s 

265 

1549 

811 

1549 

1367 

1725 

2756 

1725 

5033 

1926 



Figure 11.- Two-dimensional transient 
solid temperature distribution in char 
matrix (variable external heat flux); 

2 

where m Re .j = 68.2 kg/m -hr and 
t = 20 seconds. 
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Figure 13.- Mass flux ratio normal to 

the surface for char matrix (variable 

external heat flux); where T , # 
n 2 solid 

T fluid’ ™Ref = 68, 2 k g/ m ~ hr > ^ 
t = 20 seconds. 


Figure 12.-- Two-dimensional transient 
fluid temperature distribution in char 
matrix (variable external heat flux); 

2 

where m Ref = 68.2 kg/m -hr and 
t = 20 seconds. 


Figures 15, 16, and 17 provide the dimensionless temperature, mass flux, and 
pressure distribution, respectively, for the case of combined variable heat flux and pres- 
sure distributions together with a variable mass injection distribution at the rear bound- 
ary. The temperature distribution is very similar to that of figure 11. Once again, a 
pronounced two-dimensional in-depth flow was present; the mass flux ratio along a con- 
stant y/L line increased as the front surface was approached for y/L >0.6 and 

/ y y 

x/L < 0. 6. This, again, resulted from higher in-depth pressures associated with the 
stagnation region. 
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20 X 10" 3 



Figure 14. - Internal pressure distribu- 
tion of char matrix (variable external 

heat flux): where T , * T-. ... 

n solid fluid’ 

2 

m Re f = 68.2 kg/m -hr, and 
t = 20 seconds. 



Figure 15. - Two-dimensional transient 
temperature distribution for char 
matrix assuming thermal equilibrium 
(variable heat flux, pressure, and mass 

2 

injection); where m Ref = 67 kg/m -hr, 

Pgtag = 2. 03 x 10 3 N/m 2 , and 
t = 20 seconds. 
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Figure 16. - Internal mass flux ratio 
normal to surface for char matrix 
(variable heat flux, pressure, and mass 

2 

injection); where ihj^ = 67 kg/m -hr, 

P Stag = 2 - 03 x 1q3 N/m 2 , and 
t = 20 seconds. 


18 x 10" 3 



Figure 17.- Internal pressure distribu- 
tion of char matrix (variable heat 
flux, pressure, and mass injection); 

2 

where m R , = 67 kg/m -hr, 

1 X 61 n n 

Pgtag = 2.03 x 10 N/m , and 
t = 20 seconds. 
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CONCLUDING REMARKS 


The flow of a compressible gas through a porous matrix externally heated at the 
fluid exit surface has been investigated. The physical model was formulated to simu- 
late the flow of a pyrolysis gas formed during the degradation of an ablative material. 
The formulation of a numerical solution to the coupled energy, mass, and momentum 
relations describing this multidimensional flow of gas through a high-temperature 
porous matrix has been achieved. Two-dimensional coupled solutions were obtained 
using the strongly implicit procedure. 

In this study, the strongly implicit procedure has been successfully applied to 
the solution of a system of coupled, nonlinear partial differential equations. The 
numerical results indicate the significance of multidimensional flow for both isothermal 
as. well as nonisothermal effects in a porous material subjected to nonuniform boundary 
conditions. 


Lyndon B. Johnson Space Center 

National Aeronautics and Space Administration 
Houston, Texas, January 11, 1974 
986-15-31-04-72 
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APPENDIX A 

DERIVATION OF THEORETICAL RELATIONS 


The basic conservation relations for mass, energy, and momentum are applied 
to a compressible gas flowing through a high-temperature porous matrix. A continuum 
is assumed to exist in the porous matrix. 


CONSERVATION OF MASS 


The continuity equation for the unsteady flow of a single- species gas with a source 
term is 



(Al) 


where the superficial fluid velocity = <pv and v is velocity in a pore. 


CONSERVATI ON OF MOMENTUM . 


The conservation of momentum based on Newton's second law of motion is 


3(p f v) 

-% J +v-(p f vv ]= -vP - v-T + P f g 


(A2) 


where t is the viscous stress tensor and g is the gravitational constant. Direct 
application of this equation to flow in porous media has been discussed by Whitaker 
(ref. 12). The inability to specify exactly the flow geometry of the porous matrix and 
boundary conditions makes the solution of equation (A2) unfeasible. 

The viscous force term v • T was evaluated experimentally by Darcy (ref. 13) 
for slow, steady incompressible flow. These assumptions yield 


vP = v-t 


(A3) 
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Darcy’s results indicated that 


V-T 



(A4) 


where K is the proportionality factor. Substituting equation (A4) into the constant- 
density flow form of equation (A2) yields 


P f d? _ 

T dT"' vP " 


t 

K 


+ P f g 


(A5) 


Equation (A5) is known as the generalized Darcy Law (ref. 14). This equation 
contains both a viscous and an inertial contribution to the fluid motion. The substantial 

derivative of V contains both the temporal and the convective acceleration of the fluid. 
The basic assumption in this generalization of Darcy’s Law is the use of equation (A4) 
as representative of the viscous stress tensor. 

The validity of Darcy’s Law for slow, viscous flow has been established. How- 
ever, the use of this law at high-velocity flows does not correlate with experimental 
evidence. Equation (A6) correlates the experimental results of pressure drop and flow 
rates in a porous medium and has provided the necessary analytical correlation for 
predicting performance at other conditions (refs. 15 and 16). 


- V P 


(an + i3p f |v|)v 


(A6) 


The term /3Pj|v|v in equation (A6) is basically comparable to the operation p^V ■ v V 

in equation (A5). This modified Darcy equation (eq. (A6)) is assumed valid for use as 
the conservation of momentum relation. 

CONSERVATION OF ENERGY 


The development of the energy equation follows the approach of Bird et al. (ref. 17) 
and Clark (ref. 18). The fluid energy equation can be written, neglecting radiation and 
gravity effects, as 
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where 0 is the total internal energy. Equation (A8) is obtained by performing the in- 
dicated differentiations on the left side of equation (A7), applying the definition of the 
substantial derivative, and using the mechanical energy equation. 

~ h 

Jr - vV)] - v(V T f) - Kvv) + ^(T b - T ( ) + q'" ♦ 4 

(A8) 


where $ is viscous dissipation. Using the continuity equation (eq. (Al)) and the defini- 
tion of enthalpy, equation (A8) may be written, neglecting viscous dissipation for a 
perfect gas, as 






(A9) 


which is equation (2) in the text. 


ENERGY EQUATION FOR SOLID 


Application of the law of energy conservation to the solid state yields 


s 


3T. 

"s C p s Tf 


= v- (V T s) - Trr??( T s ' T f) + q ’” 


which is equation (3) in the text. 


REDUCTION OF GOVERNING EQUATIONS 


(A10) 


The pressure, flow rate, and temperature distributions in a porous medium are 
obtained by a simultaneous solution of equations (Al), (A6), (A9), and (A10), in conjunc- 
tion with an equation of state and the boundary and initial conditions. However, it is 
possible to reduce the number of equations by combining Darcy's Law and the continuity 
equation. 
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Assuming that the mass flux through the porous matrix can be regarded as con- 
stant for an iteration, equation (A6) can be solved explicitly for V. 


y = 

+ /3p|v|) 


(All) 


Now, using the perfect gas law 



(A12) 


And substituting equations (All) and (A12) into equation (A6) yields 


0 - 


( p \ 

PyP 

( RT f ) v * 

RT f (np + /3p|v|) 


= m' 


(A13) 


Noting that 


3P _ 1 3P 2 

at 2 P at 


(A14) 


and 


v(PvP) 



(A15) 


yields 


<t> a/p 2 \ 

i 

<1 

to 

1 

p atlT f y ' v ’ 

1 

T f (ap +j3p|v|)J 


= 2Rm"’ 


(A16) 


This is the Darcy continuity equation (eq. (1)). 
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APPENDIX B 

FINITE DIFFERENCE RELATIONS USING THE STRONGLY 
IMPLICIT TECHNIQUE 


The fluid energy equation (eq. (2)) is put in finite difference form as an illustra- 
tion of its application to the strongly implicit procedure (SIP) algorithm. The same 
procedure was followed for the Darcy continuity and solid energy equations. 

INTERIOR NODES 


Assuming for illustration purposes that m ,n = 0, q ,M = 0, and dP/dt = 0, equa- 
tion (2) can be written 


'3T, 


CP. 


Pf f\ at 


+ u 


3T f 

3x~ 



(Bl) 


Using central differences and spatially averaged thermal conductivities, equation (Bl) 
can be written as 


TI - T f /TI - TI 

i.3 i.3 _ u I i+l,j i-l,j 


At 


i ,r 


2 Ax 


+ v. . 
1.3 


'TI - TI \ 
i,j+l i,j-l _ 
^ 2 Ay / 


AX. , . + AX. . 
1+1.3 1,3 


(Ax)' 


T* 

f. 


i+1,3 


AX. - . + AX. . AX. . + AX. - . AY. . .. + AY. . 

1+1.3 1,3 , 1.3 IzlA | 1.3+1 1.3 

2(Ax) 2 2(Ax) 2 2(Ay) 2 


AY. ... + AY. . 
1.3-1 i j 3 

2(Ay) 2 


T’ + 

J M L 


AX. - . + AX. . 
1-1.3 i,3 


2(Ax)* 


T* 

i 


i-l»j 


AY. . 1 + AY. . 
1.3+1 i,3 

2(Ay) 2 


AY. . + AY. . . 
i>3 1,3-1 

lf i, i+ rr 2 (a ? ) 2 


T» 

f. . 
1 . 3-1 


v. 


1 T’ + ■ 


< pP f C x f. . 0p„ C s. . 

Pj 1.3 \ Pj 1.3 


(B2) 
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where 


AX = 


IX 

P f C 
f p» 


and AY = 


fy 

p,C 
f p f 


Formulating equation (B2) as 


AT.' . . + BT* . + CT.’ . + DT’ . + ET* , - Q. . 
i, ]-l l-l,] 1,1 1+1,3 1,3+1 ^>J 


(which is equation (16)) requires that 


v. . AY. . + AY. . , 

A - x >3 , i,3 i,3-l 

1,1 24 y 2(iy) 2 


u. AX. , . + AX. . 
B =__M + IzLl Lil 

1,1 2ix 2(ix) 2 


c t.r 


1 

At 2(ix) 2 2(Ax) 2 


AY. . + AY. . , AY. . + AY. . , v. 
1,3 1,3+1 . 1,3 1,3-1 . 1 


2(Ay)^ 


2(Ay) 


2^“ + ?p7C 


£ v 
f i P f. 


AX. , . + AX-. . u. . 
n 1+1,3 1,3 , M 

i>j 2(Ax) 2 2AX 


E. . = 


AY. - + AY, . . v, , 
J i , 3 +1 i? j 


2(Ay)' 


2Ay 


T„ / h 
f. . / v. . 

hi _ / 


,j At \0PfC s , 
l P f. 1} b 

i 


(B3) 


(B4) 


(B5a) 


(B5b) 


(B5c) 


(B5d) 


(B5e) 


(B5f) 
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BOUNDARY CONDITIONS 


As an illustration of the boundary condition requirements, consider the heated 
surface at x = 0, y = 0 (fig. 4). The SIP boundary restrictions require A = B = 0. 
Therefore, equation (B4) becomes 


CT.' . + DT.’ 1 . + ET.» . , = Q. . 
i,3 i+l,3 1,3+1 ^,3 


(B6) 


where the boundary conditions 


3T f 

W 


= 0 


or 


T* . , = T. f . , 
1,3-1 1,3+1 


(B7) 


and 


3 2 T„ 


3x 


= 0 


(B8) 


have been used. 
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